State of the Ecosystem: Indicator Processing & Visualization

Processing

d <- read.csv(file.path(raw.dir,"EcoMon Nutrient Data Through June 2018.csv"), stringsAsFactors = FALSE)

#Create data frame for mapping units to variable names
mapping <- data.frame(Units = as.character(d[1,]),
                      Var = as.character(names(d)))
mapping[mapping$Units == "" | mapping$Units == "NA",]$Units <- NA

#remove row with units
d <- slice(d,-1)

d1 <- d %>% 
  mutate(Time = Date_UTC) %>% #create Time variable
  dplyr::select(-Date_UTC,-Time_UTC) %>% #remove time, date
  gather(., Var, Value, -Latitude, -Longitude, -Time, -Depth_sampling, -Depth_station) %>% #turn wide to long while retaining lat/lon
  filter(!is.na(Value)) %>% #remove NA
  left_join(., mapping, by = c("Var")) %>% #join units 
  mutate(Longitude = as.numeric(Longitude),
         Latitude = as.numeric(Latitude),
         Time = mdy(Time)) %>% 
  filter(Latitude > 32, Latitude<50)


#Sanity check
# t1 <- d1[d1$Var == "CTDOXYMOL" ,]$Value
# t <-  d %>% slice(.,-1)
# t <- as.character(t$CTDOXYMOL)
# all(t == t1)

#Read in EPU shapefile
epu <- readOGR(file.path(gis.dir, "EPU_Extended.shp"), verbose = F) 
epu <- as(epu, "sf") #convert to sf object

if(spatial_processing){

  #Test maps
  #All EPUs
  #ggplot() + geom_sf(data = epu)
  
  #Scotian shelf
  # ss <- epu %>% filter(EPU == "SS")
  # ggplot() + geom_sf(data = ss)
  
  #get latitude and longitude for creating SpatialPointsDataFrame
  lat <-  as.numeric(d$Latitude)
  lon <- as.numeric(d$Longitude)
  coords <- data.frame(lon = lon,lat = lat)
  
  #create spdf
  spdf <- SpatialPointsDataFrame(coords = coords, data = coords,
                                 proj4string = CRS(crs))
  #convert to sf
  coords_sf <- st_as_sf(spdf) 
  
  #get intersection for mapping EPUs back to nutrient data
  epu_intersect <- st_intersection(epu, coords_sf)
  #plot(epu_intersect[epu_intersect$EPU == "MAB",])
  
  #Map back to nutrient data frame
  epu_df <- data.frame(Longitude = epu_intersect$lon,
                       Latitude = epu_intersect$lat,
                       EPU = epu_intersect$EPU)
  #join
  NE_LME_nutrients <- d1 %>% 
    left_join(.,epu_df, by = c("Latitude","Longitude"))
  
  #Select data for plotting 
  Nitr <- NE_LME_nutrients %>% filter(Var == "NITRIT.NITRAT")
  
  #Back to SOE format and specify bottom, mid-water, or surface sampling
  NE_LME_nutrients <- NE_LME_nutrients %>%
    dplyr::select(-Latitude, -Longitude) %>% 
    mutate(Value = as.numeric(Value),
           Depth_station = as.numeric(Depth_station),
           Depth_sampling = as.numeric(Depth_sampling)) %>% 
    mutate(bot_dif = Depth_station-Depth_sampling) %>% 
    mutate(surf_bot = ifelse(bot_dif <= 10, "bottom",
                           ifelse(bot_dif > 10 & Depth_sampling <= 5, "surface", "mid-water"))) %>% 
    filter(Value > 0, !is.na(EPU), !Var %in% c("BTLNBR","CASTNO","Depth_sampling",
                                             "Depth_station","STNNBR")) %>% 
    mutate(Var = paste(Var, surf_bot)) %>% 
    dplyr::select(Time, Var, Value, Units, EPU) %>% 
    group_by(EPU, Time = year(Time), Var, Units) %>% 
    dplyr::summarise(Value = mean(Value, na.rm = TRUE)) %>% 
    as.data.frame()
  
  if (save_clean){
    save(NE_LME_nutrients,file = file.path(clean.dir, "EcoMon_nutrients.Rdata"))
  }
} else {
  load(file.path(sample.dir,"sample_nutrients.Rdata"))
  load(file.path(clean.dir,"EcoMon_nutrients.Rdata"))
}

QA

Surface Nutrients

Bottom Nutrients

Processing

d <- read.csv(file.path(raw.dir, "NAO_index_1950-Oct2018.csv")) %>% slice(.,-826:-825)

#Three month running average
nao_cpc <- d %>% 
  mutate(Time = paste0(YEAR,":Q",rep(1:4, each = 3))) %>% 
  group_by(Time) %>% 
  dplyr::summarise(Value = mean(INDEX)) %>% 
  mutate(Var = "nao index", units = "unitless", EPU = "all",
         Time = yq(Time))

#annual average
nao_annual <- d %>% 
  group_by(YEAR) %>% 
  dplyr::summarise(Value = mean(INDEX),
                   Variance = var(INDEX)) %>% 
  mutate(Var = "nao index", units = "unitless", EPU = "all") %>% 
  dplyr::rename(Time = YEAR)

NAO index and index variance (Annual)

MAB Corr. Matrix

GB Corr. Matrix

GOM Corr. Matrix

Methods

To create correlation matrices, we first developed a function (shown below) to normalize and assess each time series for stationarity using the Augmented Dickey Fuller (ADF) test. If a time series were found to be non-stationary, we applied a first-order differencing step before assessing for stationarity again. If upon the second ADF test the series was still shown to be non-stationarity, it was not included in further correlation analyses rather than applying higher order differencing. This step was taken to prevent over-differencing of series. To address bias introduced by small samples sizes, only time series with N > 30 were considered for correlation analyses. Group structure in correlation matrices was defined by the magnitude of the scaled and centered first principal component.

#Normalization and differencing function

get_dat <- function(field){
  #Split out data
  time <- env[env$Var == field,]$Time
  end <- max(time)
  time <- time[1:which(time == end)]
  Value <- env[env$Var == field,]$Value
  
  if (all(is.na(as.numeric(Value))) | sd(Value, na.rm = T) == 0){
    Value <- NA #Assign as NA if not including
  } else {
    
    Value <- Value[1:length(time)]
    Value <- (Value-mean(Value, na.rm = TRUE))/sd(Value, na.rm = TRUE) #normalize
    
    #interpolate missing values
    if (any(is.na(Value))){
      Value <- approx(time, Value, n = length(Value))$y
    }
    
    #test for stationarity with Augmented Dickey Fuller
    adf <- suppressWarnings(adf.test(Value)$p.value)
    
    if (adf > 0.05){ #if non-stationary take first difference and use residuals
      mod <- arima(Value, order = c(0,1,0))
      Value <- resid(mod)
      
      adf2 <- suppressWarnings(adf.test(Value)$p.value) #check again for stationarity
      if (adf2 > 0.05){
        Value <- NA #if still non-stationary, make NA for exclusion
      }
      
    }
    
  }
  out <- data.frame(var = field,
                    value = as.numeric(Value),
                    time = time)
  return(out)
}